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Drude Weight of the Two-Dimensional Hubbard Model 
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The Drude weight of the Hubbard model on the two-dimensional square lattice is studied by 
the exact diagonalizations applied to clusters up to 20 sites. We carefully examine finite-size 
effects by consideration of the appropriate shapes of clusters and the appropriate boundary 
condition beyond the limitation of employing only the simple periodic boundary condition. 
We successfully capture the behavior of the Drude weight that is proportional to the squared 
hole doping concentration. Our present result gives a consistent understanding of the transition 
between the Mott insulator and doped metals. We also find, in the frequency dependence of the 
optical conductivity, that the mid-gap incoherent part emerges more quickly than the coherent 
part and rather insensitive to the doping concentration in accordance with the scaling of the 
Drude weight. 

KEYWORDS: Hubbard model, Lanczos method, metal-insulator transition, Drude weight, optical conduc- 
tivity 
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1. Introduction 

Metal-insulator transition is a fundamental issue in the 
condensed-matter physics. Particular interest has been 
focused on properties of metals near the transition be- 
tween the Mott insulator and the metallic phase. Many 
studies have already been made from both theoretical 
and experimental points of view. 1 

One of the quantities by which one can distinguish be- 
tween metallic or insulating phases is the Drude weight 
D, namely the coherent component of the optical conduc- 
tivity (j(uo). The system is metallic or insulating depend- 
ing on D 7^ or D = 0, respectively. In the Mott insu- 
lator, the Drude weight therefore always vanishes. If one 
dopes the Mott insulator with carriers, it becomes metal- 
lic and the Drude weight increases with further doping. 

The Hubbard model is by far the simplest microscopic 
model which can show the Mott-insulator-metal transi- 
tion. It is the most appropriate model for intensive the- 
oretical studies. Precise determination of its properties 
by some methods with no approximations contributes 
much to the fundamental understanding of the metal- 
insulator transition. In spite of the simplicity, various as- 
pects of typical many-body problems are involved in this 
model due to the presence of the Coulomb interaction. 
Many of them are still controversial even now. Among 
them is the hole-doping dependence of the Drude weight 
near the transition between the Mott insulator and the 
metal. Primary interest of this paper is the dependence 
of the two-dimensional Hubbard model for small doping 
concentration; this is because the Mott insulator of this 
model at half filling is considered to become metallic in 
a nontrivial way when the doping is switched on and be- 
cause this model has been extensively studied as a model 
of high-T c cuprates. 

Generally, it is not so easy to calculate such a dynam- 
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ical quantity as the Drude weight of the strongly corre- 
lated system without approximations. One possible way 
is to apply the exact diagonalization method to finite- 
size clusters. The first work along this line concerning 
with the two-dimensional Hubbard model was reported 
by Dagotto et al? Based on the treatment on the clusters 
up to 16 sites, they successfully captured a rough feature 
of the Drude weight as a function of the hole doping con- 
centration 5. They tentatively concluded that D oc J, 
although the Drude weight sometimes becomes nega- 
tive with large absolute values especially near the metal- 
insulator transition point. Such negative Drude weights 
occur as a consequence of the finite-size effect. In order 
to get rid of such a finite-size effect, Nakano and Imada 
examined the properties under various boundary condi- 
tions in the system of the same size with 16 sites. 3 They 
consequently found the concave behavior of the Drude 
weight as a function of S in the small-c) region, which 
implies the dependence D oc S 2 . Recently, Tohyama and 
his co-workers treated larger clusters of 18 and 20 sites 
under the periodic boundary condition. 4 Their result of 
the Drude weight seems to suggest D oc S again. Thus, 
the issue of the (^-dependence of D is still controversial 
at present. 

Under these circumstances, we examine the finite-size 
effects of the Drude weight of clusters with various shapes 
and sizes up to 20 sites under various boundary condi- 
tions very carefully and do our best in minimizing the 
finite-size effects when we calculate the Drude weight of 
the finite-size clusters. The purpose of the study in this 
paper is to capture thereby a reliable exponent for the 
^-dependence of the Drude weight near the transition 
between the metal and the Mott insulator. Our present 
results suggest that D oc S 2 . 

This paper is organized as follows. In the next section, 
the model Hamiltonian and the method of numerical cal- 
culation are introduced. Section 3 is devoted to the ex- 
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amination of the finite-size effect in the non-interacting 
case. In Section 4, results in the interacting case are re- 
ported and discussed. Summary and conclusion are given 
in the final section. 

2. Model and Method 

We consider the two-dimensional Hubbard model. Its 
Hamiltonian is given by 



H 

7~^hop 
T^int 



7Yh 



Hh 



U 



op 



(1) 



where the creation (annihilation) of an electron at site 
i with spin a is denoted by c\ G (q j0 -) and is the 

corresponding number operator, namely = c\ a Ci j(T . 
In this paper, we treat only the nearest-neighbor hop- 
ping on the two-dimensional square lattice. Shapes for 
finite-size clusters will be illustrated later. Energies are 
measured in units of t\ therefore, we set t = 1 hereafter. 
The system consists of N e electrons on 7V S atomic sites. 
The hole doping concentration is a controllable param- 
eter given by 5 = 1 — n, where the electron density is 
defined as n = N e /N s . 

We perform numerical diagonalization of finite-size 
clusters of the model (1) using Lanczos algorithm. This 
method is non-biased against the effects of Coulomb 
interaction. The method is also useful to obtain the 
ground-state wave function with the ground-state energy 
E g very precisely. On the other hand, the disadvantage 
of this method is that available system sizes are lim- 
ited to be small. This is because the dimension of the 
Hilbert space grows exponentially with respect to the 
system size. 5 In order to overcome this disadvantage, we 
perform parallel calculations in the computer systems 
of the distributed- memory type. Thereby we can treat 
arbitrary shapes of clusters up to N s = 20 if we use ap- 
propriate supercomputer systems. We also employ the 
continued-fraction expansion method, 6 which provides us 
with dynamical quantities including the optical conduc- 
tivity and the Drude weight. The detailed procedure will 
be explained after the examinations of cluster shapes and 
boundary conditions. 

3. Finite-Size Effect for U = and Boundary 
Condition 

In this section, we examine how finite-size effects oc- 
cur in various types of boundary conditions of clusters 
for a given N s . When one considers the system of two- 
dimensional square lattice by using a finite-size cluster, 
regular squares under the periodic boundary condition 
are usually adopted as in refs. 2 and 4. Such regular 
squares for N s = 16, 18 and 20 are depicted in Fig. 1 (a), 
(b) and (c), respectively. The merit of adopting these 
squares is that the x and y directions are equivalent. 
When one treats such a quantity as the Drude weight, 
however, calculations based on the regular squares under 
the periodic boundary condition reveal a large finite-size 




Fig. 1. Finite size clusters for N s = 16, 18 and 20. The blue 
and green vectors represent two translational vectors for each 
cluster. The atomic sites are denoted by the red circles inside 
the quadrangle formed by these two vectors and the other edges 
denoted by black solid lines. 



deviation from the value in the thermodynamic limit. 
In ref. 2, a remarkable finite-size effect was reported. In 
view of this situation, effects of other types of boundary 
conditions, namely the antiperiodic boundary condition 
and the mixed boundary condition, were studied in ref. 
3 to see whether one can reduce the finite-size effects 
for the fixed N s or not. Consequently, the wrong nega- 
tive large Drude weights were successfully resolved. Note 
here that when one imposes the mixed boundary condi- 
tion, the two lattice directions are no longer equivalent 
and the average of the two directions gives a reasonable 
result. This suggests that adhering to the equivalence of 
the two directions is not necessary. The two directions 
become inequivalent when one relaxes the condition that 
these quadrangles of the clusters have the right angle at 
their vertices. Therefore we rather additionally examine 
the parallelogram clusters depicted in Fig. 1 (d) and (e) 
for N s = 16 and 18 in place of squared-shaped clusters, 
respectively. An important point is that a set of possible 
wave number vectors corresponding to the clusters (d) 
and (e) is different from the one corresponding to the 
clusters (a) and (b), respectively. The system is never- 
theless still bipartite. For N s = 20, on the other hand, 
it is impossible to preserve the bipartite nature when 
one distorts the cluster (c) so that it realizes a parallel- 
ogram cluster with a different set of possible wave num- 
ber vectors. Then we study only one type of cluster for 
N s = 20. Therefore we consider these clusters (a)-(e) up 
to N s = 20. 

Next, we show how to select the appropriate cluster 
under the appropriate boundary condition for a given 
N s . For this purpose, let us first see the ground-state en- 
ergy E g of the clusters (a)-(e) at half filling for [7 = 0, 
presented in Table I. Note that this quantity is propor- 
tional to the Drude weight for U = 0. A noticeable point 
of the results of clusters (a), (b) and (c) under the pe- 
riodic boundary condition is that the relative difference 
(i^oo — E g )/Eoo does not monotonically converge as N s 
is increased, where Eqq denotes the ground-state energy 
in the limit of N s — > oo. This suggests that the system 
size 7V S up to 20 is not free from the finite-size effect as 
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far as we stick to the periodic boundary condition. We 
have thus to be careful in choosing an appropriate clus- 
ter under an appropriate boundary condition for a given 
N s . Let us now compare the difference (E^ — E g )/E OQ of 
all the clusters (a)-(e) with various boundary conditions. 
One can find the best cluster and boundary condition 
from among the cases shown in Table I so that the abso- 
lute value of the difference is the smallest. For N s = 16 
and 20, two boundary conditions give the smallest one. 
Sets of possible wave number vectors of the two cases 
for each N s are equivalent. Thus, the same result will be 
obtained if one selects either of the two. 

Let us compare the ^-dependence of the non- 
interacting E g between two cases, the one of the most 
appropriate cluster under the most appropriate bound- 
ary condition with the best E g at half filling and the 
other of the regular square under the periodic boundary 
condition; results are depicted in Fig. 2. One can easily 
observe the results of the former case fall closer to the 
curve of the thermodynamic limit in the whole region of 
5, irrespective of the size N a . On the other hand, the lat- 
ter regular-square case under the periodic boundary con- 
dition gives significant overestimates particularly at (N s , 
N e ) = (18, 18) and (20, 18). This strongly suggests that 
the Drude weights of these cases will be possibly overesti- 
mated even when U is considerably large. It also suggests 
that if the intrinsic Drude weight shows the concave be- 
havior as a function of S for large Z7, this overestimation 
will possibly cancel this behavior. Figure 2 also provides 
another important finite-size effect that finite-size data 



Table I. List of ground-state energy per site in the non- 
interacting case U = at half filling for various shapes of clusters 
and various boundary conditions together with relative difference 
from the energy in the thermodynamic limit. Head BC1 (BC2) 
denotes the boundary condition along the green (blue) transla- 
tional vector in Fig. 1. P and A represent the periodic boundary 
condition and the antiperiodic boundary condition, respectively. 
Note here that the boundary conditions giving the same energy 
for each cluster correspond to the identical set of possible wave 
number vectors. Marked cases in last column with head PW show 
the selected clusters under the selected boundary condition in 
which we actually perform calculations in the present work. 



cluster 


BC1 


BC2 


-E g /N s 


(Eoo — Eg)/ Eoo 


PW 


(a) 


P 


P 


1.5 


-0.0747 






A 


P 


1.70711 


+0.0530 






P 


A 


1.70711 


+0.0530 






A 


A 


1.41421 


-0.1276 




(b) 


P 


P 


1.77778 


+0.0966 






A 
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1.53960 


-0.0503 






P 


A 


1.53960 


-0.0503 






A 


A 


1.33333 


-0.1775 




(c) 


P 


P 


1.69443 


+0.0452 






A 
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1.65065 


+0.0182 


• 
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1.65065 


+0.0182 
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1.52169 


-0.0613 




(d) 
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1.70711 


+0.0530 






A 
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1.63099 


+0.0061 


• 
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1.5 


-0.0747 






A 


A 


1.63099 


+0.0061 




(e) 


P 
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1.72417 


+0.0636 
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1.64519 


+0.0148 


• 
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1.51621 


-0.0647 
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1.44675 


-0.1076 
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Fig. 2. Ground-state energies as function of 5 in non-interacting 
case for given N s presented in the form of — Eg/ (4iV s ). This quan- 
tity corresponds to D for U = 0. Green triangles represent data 
for the regular square under the periodic boundary condition. 
Red circles represent data for the appropriate cluster under the 
appropriate boundary condition based on the analysis in Table 
I. The solid curve denotes the corresponding quantity for the 
thermodynamic limit. 



show nonzero curvature in their ^-dependence only at 5 
corresponding to the closed electronic shell structure. At 
S for the open electronic shell structure, on the other 
hand, E g follows linear dependence. The linear depen- 
dence from this origin is likely to survive more or less 
even for U > 0. This is, in fact, confirmed in the quantum 
Monte Carlo study at nonzero U. 7 Thus, it is difficult to 
judge whether the linear ^-dependence revealed from nu- 
merical studies in the interacting case in the region of the 
open shell structure is intrinsic or is simply attributed to 
this particular electronic structure. In view of this situ- 
ation, it is desirable to confine our analysis to the data 



4 J. Phys. Soc. Jpn. 



Full Paper 



H. Nakano, Y. Takahashi and M. Imada 



at S for the closed shell structure in order to extract the 
intrinsic behavior. In this paper, therefore, we take the 
appropriate clusters under the the appropriate boundary 
condition for each N s (shown in Table I) and adopt data 
for the filling forming the closed shell structure. 

In the selected clusters and boundary conditions, the 
equivalence of the two directions denoted by a (= x,y) 
does not hold, as we have already mentioned in the above. 
Then, we take the direction average according to the pro- 
cedure used in ref. 3. The optical conductivity is given 
by a(uj) = 2ttDS(uj) + a reg (cc;). Here, the second term is 
the regular part of the optical conductivity, which is also 
called the incoherent part. The regular part a reg (cc;) is 
obtained by the average a reg (cc;) = [a r x eg (u) + a r y eg (u)}/2; 
each of <7 Y ^ g {uo) given by 

mm? 



reg 



- v 



-S(lo -E e + Eo), (2) 



is calculated by the continued-fraction expansion 
method. Here, j a is a current operator along the a- 
direction defined as 



3a — i'5~*' J t( c i,er C i+S a ,cr c\ + 5 ata Ci,a) , 
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Fig. 3. Ground-state energies as function of S in interacting cases 
presented in the form of — Eg/(4N S ). The solid curve denotes 
— Eg/(4N S ) for U = for the thermodynamic limit. Closed cir- 
cles and closed triangles denote the present results from exact 
diagonalization calculations for U = 4 and U = 35, respec- 
tively. Open symbols (triangles, squares, reversed triangles and 
diamonds) denote results of TVs = 6 x 6, 8 x 8, 10 x 10 and 
12 x 12 for U = 4 from the quantum Monte Carlo simulations, 
respectively. 7 



where 8 a is the unit vector along the a-direction, \£) rep- 
resents an eigenstate with the energy eigenvalue Eg. Note 
that the ground state is described by i = 0. One can ob- 
tain the Drude weight D from the combination of (t(uj) 
and the sum rule 



f 

Jo 



(4) 



where — AK represents the kinetic energy per site, namely 
— 4K = (7Yhop)/^Vs- For U = 0, the incoherent part is 
absent; therefore one has D = K. It is widely known that 
for U > 0, the incoherent part appears. When U is large 
enough, one can easily recognize in a(uo) the responses 
toward the upper-Hubbard band. In such a case, one can 
definitely determine a frequency just below the upper- 
Hubbard band to be cj c ; we also calculate definite integral 



eff 



__ 1 r 

^ Jo 



(5) 



which is called an effective carrier density. 

As for the validity of extending the present procedure 
to cases with finite U, there is no convincing theoretical 
argument on the condition of reducing the finite-size ef- 
fect. Though not reliable in a strict and rigorous sense, it 
is at least better to start from situations where the effects 
are reduced in the case of U = from the continuity of 
behavior as the function of U. Another evidence for the 
validity of the present procedure is given from our exam- 
ination of E g for U > 0; the results are depicted in Fig. 3. 
Although all the possible data of N s = 16, 18 and 20 are 
shown there with the same closed symbol for each case 
of U > 0, (5-dependence of E g is fairly smooth. This fact 
implies that finite-size deviations of the present data are 
small. One can actually confirm that the present results 
for U = 4 agree well with the results from the quan- 
tum Monte Carlo simulations 7 on larger systems up to 
7V S = 12 x 12. Therefore, it is reasonable to presume that 



the finite-size effects in D and TVeff based on the present 
procedure are also small if the smooth (5-dependence of 
D and N e ^ is confirmed irrespective of the size 7V S . 

In this paper, we take three cases of U = 10, 20 and 
35. In the system with a nonzero U weaker than U — 10, 
it is difficult to distinguish the upper-Hubbard band and 
the lower-Hubbard band in the responses of ct(uj); the 
estimation of 7V eff is not available. On the other hand, 
in the system for U > 35, the ground state is possibly 
ferromagnetic near half filling. Once the ferromagnetism 
appears, it is too difficult to study the scaling of D near 
half filling. 

4. Results for Interacting Case and Discussions 

We show numerical results for the (^-dependence of N e ^ 
connected by lines for interacting cases in Fig. 4. The 5- 
dependence of the Drude weight is later shown. Although 
the data consist of those of various system sizes for each 
case of U, ^-dependence of them is fairly smooth, inde- 
pendent of the value [/, as shown in Fig. 4(a). It clearly 
indicates that deviations of our numerical results due to 
the finite-size effect are very small for every U that we 
treat. One can also observe that as 6 is decreased to- 
ward the half-filling case 5 = 0, N e s goes to a very small 
value for each of non-zero U. The values of iV e ff at half 
filling are very close to the corresponding D. This is be- 
cause no significant responses within < uo < uo c are 
observed in the optical conductivity; this behavior has 
already been reported in all the previous works 2 4 and 
will be confirmed later in the present study. Then, we will 
later come back to these values of N e ^ at half filling in 
the discussion of the Drude weight. In order to examine 
the 5-dependence of N e ^ near half filling, Fig. 4(b) de- 
picts the 5-dependence of N e ^/S instead of raw N e ^. One 
can observe the linear behavior irrespective of the value 
of U in the wide region of S. They extrapolate to finite 
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Fig. 4. 5-dependence of N e ^ in (a) and N e ^/5 in (b). Closed tri- 
angles, open squares and closed circles denote data for U = 10, 
20 and 35, respectively. In (a), total weight for U = is also 
depicted by the solid curve for comparison. 




Fig. 5. 5-dependence of Drude weight for U = 10, 20 and 35. 
Notations for closed triangles, open squares and closed circles 
are the same as those in Fig. 4. For comparison, solid curve 
represents D for U = 0, which is the same as the total weight. 



intercepts of the N e ^/S axis, resulting in the following 
dependence 



TVeff oc 5, 



(6) 



near half filling 6 = 0. Note that the above dependence 
(6) is in agreement with all the previous studies. 2 4 

We depict ^-dependence of the Drude weight in Fig. 5. 
Note that the Drude weight at 6 = vanishes for U > 
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Fig. 6. Log-log plot of D and N e ff as function of S for U = 35. 
Circles and triangles denotes D and N e ^, respectively. 



in the thermodynamic limit. All the values for various U 
are so small, less than 0.003, that we cannot distinguish 
their symbols from the origin as shown in Fig. 4(a) and 
5, where our data of the half-filling case for N s = 20 are 
plotted. For N s = 16, we obtain D - 0.009 for U = 10 
at most; since its difference from the data of N s = 20 is 
only within the symbol size, it is not shown in the figure. 
On the other hand, ref. 4 reported the Drude weight of 
about 0.02, which is substantially larger than ours, for 
the cluster (c) with 7V S = 20 under the periodic boundary 
condition. Our present result is closer to the correct value 
of zero in spite of the smaller size 7V S , which indicates su- 
periority of our present calculations in reproducing well 
the insulating state at half filling. 

In the metallic region for 6 > 0, all the data obtained 
from studies on various system sizes are presented in 
Figs. 4(a) and 5 according to the discussion in section 
3. Recall here that refs. 2 and 4 reported anomalous and 
discontinuous deviations of D at some particular hole 
concentrations, resulting from a finite-size effect. 11 On 
the other hand, the present data for all the values of S 
in Fig. 5 show smooth variation on the whole, though 
they consist of mixture for various N s . It indicates that 
the finite-size effect is well eliminated in our estimates. 
As for the (5-dependence of the Drude weight, the convex 
behavior is apparent in all the region of 6 for U = 10. For 
U = 20 and 35, on the other hand, the concave behavior 
occurs near half filling. The change in this ^-dependence 
with increasing U will be discussed below. 

In order to confirm the concave behavior of D, we plot 
the data of D on the logarithmic scale in Fig. 6 together 
with 7V e ff for U = 35. Our data of D near half filling 
indicate the dependence 

D oc S 2 . (7) 

This dependence is in contrast to the linear dependence 
of 7V e ff of eq. (6). For more detailed analysis of the depen- 
dence of D in the region of small (5, let us next examine 
the (5-dependence of the ratio D/N e ^. The correspond- 
ing quantity was also studied in the two-dimensional t-J 
model. 12 The present results are shown in Fig. 7 in the 
form of (D/N e s)/5. For U = 10 and 20, as S is decreased 
toward half filling, the quantity (D/N e s)/5 shows steep 
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Fig. 7. 5-dependence of D/(N e ^5) and its inverse. Notations for 
closed triangles, open squares and closed circles are the same as 
those in Fig. 4. 



Fig. 8. 5-dependence of D/S. Notations for closed triangles, open 
squares and closed circles are the same as those in Fig. 4. Broken 
lines are the fitting lines in 0.25 < S < 0.5. Dotted lines are guides 
for eyes in the assumption of D oc S 2 . 



increase. For U = 35, on the other hand, no such behavior 
is observed. The increasing behavior for U = 10 and 20 is 
examined in the plot of the inverse of (D/N e s)/S shown 
in the inset of Fig. 7. Our data for U = 10 and 20 reveal 
linear dependence in a wide region of S. For U = 35, a 
similar linear dependence is observed in 0.2 ^ 5 ^ 0.5 
together with a deviation from it in the smaller-^ region. 
Each of the linear extrapolation of plots intersects at a 
non-zero intercept with the axis of ordinates in the limit 
of S — > 0. If one takes its weak (5-dependence for U = 35 
into account, it is natural to consider that [(D/Nq^/S] -1 
for various U stays finite in the limit of S — > as common 
behavior; there is no sign of vanishing [(D / N e ^) / S] -1 
from our data of finite size clusters. The above behav- 
ior, together with eq. (6), thus leads us to the scaling 
(7) valid not only for U = 35 but also for U = 10 and 
20. This analysis is particularly important since it sug- 
gests the validity of the scaling (7) even for U = 10 in 
spite of the fact that the ^-dependence of D for this U 
in Fig. 5 does not show concave behavior which appears 
for U = 20 and 35. 

Let us now discuss the present result from the view- 
point of the scaling theory developed by Imada, a general 
theory of the metal-insulator transition. 13 This helps us 
to justify the dependence (7) in relation to other behav- 
ior of the two-dimensional Hubbard model. According to 
the scaling theory, the Drude weight shows the following 
critical behavior near the metal-insulator transition, 

D oc S 1+ ^' d , (8) 

where z denotes the dynamical exponent and d is the 
spatial dimension. Since d = 2 in the present problem, 
one obtains z = 4 from (7). The scaling theory also pre- 
dicts the critical behavior of the compressibility near the 
transition, 

k ex ^- z ' d . (9) 

Substitution of z = 4 and d = 2 in the above dependence 
of n leads to n ex c) -1 . The diverging behavior of n was 
reported by the quantum Monte Carlo (QMC) simula- 
tions 7 ' 14 performed on the same model with larger sys- 



tem sizes not yet treated by exact-diagonalization stud- 
ies. The QMC calculation 15 also demonstrated that the 
localization length suggests the divergence^ ex \ji—ji c \~ v 
with v = 1/4 as the chemical potential approaches the 
charge gap ji c from the insulating side. The scaling the- 
ory indicates zv = 1, which also leads to the exponent 
2 = 4. Therefore, our present result is consistent with 
these QMC results in views of the scaling theory. Re- 
cent study of the metal-insulator transition 16 18 clari- 
fied that the exponent z = 4 appears at the marginal 
quantum critical point between the continuous metal- 
insulator transition with T c = and the discontinuous 
transition with nonzero T c . The appearance of the same 
marginal exponent z = 4 was also found in the one- 
dimensional Hubbard model with next-nearest-neighbor 
hopping. 10 ' 19 The unified understanding of the marginal 
point occurring in the one- and two-dimensional cases is 
the subject for future studies. Although it is difficult to 
determine the marginal quantum critical point precisely, 
the present results suggest that the critical region of the 
marginal point is wide extending to all the values of U 
we studied. This allows us to capture successfully the 
distinct signature characteristic to the marginal point. 

Next, let us consider the width of the critical region 
of S in which the behavior (7) is observed. For this 
purpose, we depict the ^-dependence of D/S in Fig. 8. 
In 0.25 ^ S ^ 0.5, one can observe [/-irrelevant lin- 
ear behavior. In Fig. 8, we fit our numerical data in 
0.25 < 5 < 0.5 with broken lines. For U = 20 and 35, 
our data near half filling digress from the broken lines. 
In particular, for U = 35, two points shows clear de- 
viation toward the origin in Fig. 8, while one point for 
U = 20. Here, let us draw the dotted guide lines for 
these digressing data by assuming the present result (7). 
The dotted and the broken lines cross with each other for 
each U. One can regard the crossing point for each [/asa 
rough estimation of the boundary of the critical region. 
Then, the estimated boundaries are S ~ 0.14 and 0.22 
for U = 20 and 35, respectively. Though crude are these 
estimations, this is the first study to report the bound- 
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ary. Calculations of systems with larger N s with proper 
treatment of finite-size effects will make the precise es- 
timations possible in future studies. For U = 10, there 
are no data points showing digressing towards the origin. 
If the S 2 dependence of D, i.e. eq. (7), is true, then the 
above fact suggests that the critical region for U = 10 is 
narrower than S < 1/9. This is consistent with the rapid 
decrease of the estimated boundary as U decreases from 
35 to 20. The doping concentration at the extrapolated 
boundary for U = 10 may be very small. 

Let us discuss the frequency dependence of a reg (cc;) 
obtained by eq. (2) for the insulating and metallic cases 
where the system is in the critical region; we present our 
results in Fig. 9. Figures 9(a) and 9(b) depict the results 
for insulating cases. The structures of a reg (uo) between 
N s = 16 and N s = 20 are similar for each U. This fact 
suggests that finite-size effects are rather irrelevant for 
the structure of <j reg (uo). For each J7, broad responses 
corresponding to the transitions to the upper-Hubbard 
band appear around uo = U. For U = 10, a sharp peak 
emerges at the lower edge of the band. According to ref. 
4, this sharp peak is ascribed to the spin-polaron forma- 
tion in the photoexcited state, and the weight of the peak 
becomes smaller for larger N s . On the other hand, our 
sharp peak for U = 10 does not show any significant N s 
dependence, which implies that it will survive in the ther- 
modynamic limit, although the peaks get rapidly weaker 
for larger U = 20 and 35. The broad responses are unaf- 
fected by the change in U between U = 20 and 35 except 
that the center of responses moves with U. Figures 9(c), 
9(d) and 9(e) depict the metallic cases; one can observe 
not only the high-frequency spectrum due to the transi- 
tions to the upper-Hubbard band but also a pronounced 
structure in the low-energy region below the Mott gap. 
First, let us see the transitions to the upper-Hubbard 
band. By switching S on, these transitions gradually lose 
their weights by transferring them to the lower-energy 
structure, as already been reported in refs. 2, 3 and 4. 
Next, let us see the structure in the low- frequency re- 
gion. We find strong responses growing rapidly towards 
uo = with decreasing uo. It is noticeable that the shape 
seems to be a common property irrespective of U and S 
within the present parameter ranges. The same shapes 
are produced in the present results of both N s = 18 and 
20. The cj-dependence of the shape reminds us of the 1/uo 
tail and the mid-IR peak observed in the experiments of 
high-T c compound. 20 

Finally, let us make comparisons with experimental re- 
sults of optical conductivities. Recall here that the coher- 
ent and incoherent parts are not definitely separated in 
the experimentally observed optical conductivity. Thus, 
it is not so easy to make direct comparison between ex- 
perimental and numerical <r reg (c<;). Theoretically, on the 
other hand, one has to introduce the broadening of delta 
functions with some widths in dealing with finite sys- 
tems as the present ones. Keeping these circumstances 
in mind, we consider the frequency dependence of a(uo) 
including the Drude part; Figure 10 depicts results of 
U = 10, typical value of the Coulomb interaction for 
most of the high-T c materials. For (5 = 0, the present 
spectrum reveals a clear increase as a; — > 0, even when 
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Fig. 9. Frequency dependence of regular part of the optical con- 
ductivity in (a) for N s = 16 and N e = 16, namely 6 = 0, (b) 
for 7V S = 20 and N e = 20, namely 6 = 0, (c) for N s = 18 and 
N e = 16, namely 6 = 1/9, (d) for N s = 18 and N e = 16, namely 
6 = 1/9 and (e) for N s = 20 and N e = 16, namely 6 = 0.2. Blue, 
green and red lines denote the cases for U = 10, U = 20 and 
U = 35, respectively. Delta functions are broadened with width 
of 77 = 0.05. 



the system is in the insulating state. This comes from 
the tiny D due to the finite-size effect. For 5 > 0, on 
the other hand, the Drude weight is considerably large 
in comparison with this tiny D at half filling; there oc- 
curs a huge peak at uo = as the coherent component. 
Its tail in Fig. 10 appears due to the broadening of width 
77 = 0.2. Note here that even the tail is sizable in compar- 
ison with the intrinsic responses of <j reg (o;) in the inner- 
gap region. This means that in order to understand op- 
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Fig. 10. Frequency dependence of optical conductivities for U = 
10. The insulating case of iV s = 20 and N e = 20 is denoted by 
the sky-blue dotted line. Red, deep-blue, green and black solid 
lines correspond to the metallic cases of TVs = 18 and N e = 16, 
TVs = 20 and N e = 16, N s = 16 and N e = 12 and N s = 18 and 
N e = 12, respectively. Delta functions are broadened with width 
of rj = 0.2. For comparison, the coherent part in <j(lj) is also 
presented for the two cases of N s = 18 by the broken lines of the 
same colors. 

tical responses of the system, it is important to capture 
both the behavior of the coherent component and that of 
the incoherent one included in the total cr(uo) simultane- 
ously. In Fig. 10, qualitatively the same change by doping 
in the transitions to the upper-Hubbard band is seen as 
that mentioned above for U = 20 and 35. The structure 
in this uo region should be associated with the charge 
transfer from the copper d to the oxygen p states in the 
case of the cuprates, hence it is, in narrow sense, dif- 
ferent from the present case of the Mott-Hubbard type. 
Aside from the character of the gap, i.e., either it is the 
Mott-Hubbard gap or the charge transfer gap, the evolu- 
tion of doping observed in the experiments are correctly 
captured in our numerical results, the spectral weight 
transfer from the higher to lower energies as well as the 
collapse of the sharp peak at the low-energy edge. 

The evolution of the low-energy structure with doping 
contains both the growth of the Drude weight and the 
development of the incoherent mid-IR structure. The to- 
tal sum of these two contributions is the effective carrier 
density iV e ff, whose weight is proportional to the dop- 
ing concentration. In the doping process of the evolu- 
tion from the insulator to low doping, the mid-IR in- 
coherent part below the Mott gap grows quickly in our 
calculated results and is already prominent even in un- 
derdoped systems (see the result at S ~ 0.111), with a 
rather insensitive dependence on the further doping to 
the overdoped systems (see S = 0.333). On the contrary 
to the mid-IR part, the Drude part grows slowly at low 
doping, while it starts growing quickly from moderate 
to over doping, which is the origin of the scaling given 
by D oc 5 2 rather than D oc S. The total sum resulting 
in iV e ff evolves smoothly and linearly with S because of 
the compensation of these two contributions. The fact 
that the Drude part remains only a small fraction of N e s 
whereas the mid-IR structure rapidly grows from zero to 



low dopings is indeed common to all the high-T c cuprates. 
It is even more universal in the transition metal oxides in 
general. Our calculated result for the further doping pro- 
cess from moderate to over doping is also consistent with 
the universal trend of the experimental results in transi- 
tion metal oxides including the high-T c cuprates, where 
the mid-IR structure is absorbed to the Drude weight. 
The S 2 scaling of the Drude weight is tightly related to 
this profound and universal trend in the doping process 
of the whole spectral weight in the optical conductivity. 

5. Conclusion and Discussion 

We have investigated the Drude weight of the Hub- 
bard model on the square lattice by means of large-scale 
parallelized exact diagonalizations. The hole-doping de- 
pendence of the Drude weight near the transition point 
between the metal and the Mott insulator has been a 
controversial problem. We have very carefully examined 
finite-size effects and we have successfully captured the 
intrinsic behavior of the Drude weight. Our results sug- 
gest that the Drude weight is scaled by D oc J 2 , which 
is consistent with results from various kinds of analy- 
sis. In particular, it is important that this scaling agrees 
with the quantum Monte Carlo results through the scal- 
ing theory of the metal-insulator transition. The behav- 
ior D oc S 2 is a characteristic feature near the marginal 
quantum critical point of the transition between the Mott 
insulator and the doped metal in two-dimensional sys- 
tems. The restructuring of the electronic structure upon 
doping into the Mott insulator first causes the spectral 
weight transfer from the weight above the gap mainly 
to a mid-IR incoherent component below the gap. The 
spectral weight is further progressively transferred from 
the mid-IR component to the Drude weight. This large- 
scale restructuring seen in the frequency and doping de- 
pendences of a(u) in our result is consistent with the 
experimental results on transition metal oxides includ- 
ing high-T c cuprates. The 5 2 scaling of the Drude weight 
comes from this underlying physics involving the high 
energy scale. 

AtU = 0, the Hubbard model we studied has a perfect 
nesting just at half filling and the Fermi level coincides 
with the van Hove singularity. The effect of deviation 
from this perfect nesting by introducing the hopping to 
further sites (such as the next-nearest-neighbor hopping 
t') is left for separate study. The QMC study 14 indicates 
that the scaling behavior of the compressibility does not 
sensitively depend on the presence or absence of t' ', which 
implies that the present scaling of the Drude weight also 
holds. Next, let us recall the marginal point with z = 4 
in the one-dimensional Hubbard model with the next- 
nearest-neighbor hopping 10 ' 19 mentioned above. In this 
case, a singularity in the density of states is absent at 
half filling at U = because the next-nearest-neighbor 
hopping is small. This fact also suggests that the metal- 
insulator transition with z = 4 occurs irrespective of such 
a singularity in the density of states. 

The finite-size effects in such a quantity as the Drude 
weight are considerably large. Thus, the careful exam- 
ination of the finite-size effects is necessary in possible 
ways as much as we can. If one calculates the quantity 
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without considering the effects, obtained results could 
severely suffer from superficial finite-size effects. By such 
results, unfortunately, the incorrect conclusion would be 
derived. In this meaning, the present problem of the 
Drude weight is an instructive example. As a way to re- 
duce the finite-size effect, the present work employs ap- 
propriately chosen clusters under the appropriate bound- 
ary conditions. Thereby, we have captured the critical 
behavior of D oc S 2 near the half- filling insulator. Av- 
eraging over twisted boundary conditions 21 may be an- 
other possibility for the reduction. Quite recently such 
an investigation has been done in the N s = 20 cluster of 
the t-J model plus three-site terms. 22 Although this work 
gives a result of the (5-dependence of D, it is so subtle that 
one cannot conclude either of D oc 5 2 or D oc 5. In the 
not-too-distant future, computers will be even more pow- 
erful; numerical diagonalizations of the Hubbard model 
with larger system sizes would be possible. Such calcula- 
tions could clarify the validity of the present procedure 
for reducing the finite-size effects in D and provide us 
with a more definite result for the critical behavior of D. 
We hope that the knowledge in the present study will 
serve for such future studies. 
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